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PREWHITENING  OF  COLORED  NOISE  FIELDS  FOR 
DETECTION  OF  THRESHOLD  SOURCES 

INTRODUCTION 


High  resolution  direction  finding  methods  such  as  Multiple  Signal  Classification 
(MUSIC),  1  Enhanced  Minimum  Variance  Distortionless  Response  (EMVDR),^  and 
Minimum  Norm  (Min-Norm)^  compute  angular  response  as  a  function  of  cosine  of  bearing 
angle.  These  methods  associate  peaks  in  the  response  function  with  source  locations, 
providing  direction  of  arrival  estimates.  A  response  peak  that  does  not  correspond  to  an 
existing  source  indicates  a  false  alarm.  Rooting  techniques  for  direction  finding  such  as 
ROOT-MUSIC,4  Estimation  of  Signal  Parameters  Via  Rotational  Invariance  Techniques 
(ESPRIT), 5  and  the  Pisarenko  Harmonic  Decomposition  (PHD)6  provide  source  locations 
directly,  without  computation  of  a  response  function,  by  constructing  a  polynomial  whose 
roots  correspond  to  arrival  angles.  Any  angle  that  does  not  correspond  to  an  existing  source 
indicates  a  false  alarm.  Typically,  these  algorithms  provide  a  large  set  of  candidate  arrival 
angles  that  includes  true  source  locations  and  spurious  ones.  Although  methods  exist  for 
determining  the  true  source  locations  from  the  candidate  scij  they  require  an  accurate  set  of 
candidate  angles  with  as  few  false  alarms  as  possible. 

Both  types  of  direction  finders  exhibit  high  false  alarm  rates  in  the  presence  of  spatially 
colored  ambient  noise.  In  their  traditional  model  formulation,  direction  finders  assume  that 
only  spatially  white  noise  exists.  When  an  extended  or  continuous  colored  noise  is  present 
in  the  acoustic  field,  low  SNR  or  threshold  sources  may  be  obscured  from  detection.  Also, 
the  estimate  of  a  source  location  for  a  source  that  exists  within  the  spatial  bandwidth  of  the 
colored  noise  may  be  biased.  Thus,  in  the  presence  of  known  nonwhite  noise,  high  false 
alarm  performance  necessitates  the  use  of  prewhitening  algorithms.  8  Direction  finders  also 
require  a  reliable  estimate  of  the  number  of  sources  present  in  the  observed  noise  field. 8 
Because  nonwhite  noise  skews  or  inflates  this  estimate,  prewhitening  is  required. 
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MOTIVATION  FOR  PREWHITENING 


Define  x  as  the  N-dimensional  sensor  output  data  vector  from  an  N-element  uniform  line 
array.  Denote  the  covariance  matrix  of  x  by 

R  =  E[xxh]  ,  (1) 

where 

*T  *  [*!  *2  -  x«] .  O 

E[]  represents  the  expectation  operator,  and  H  denotes  complex  conjugate  transpose.  Model  R 
with  the  following  structure^ 

R  =  o^P  +  o^Q  +  <#N  •  (3) 

where  P  corresponds  to  the  discrete  source  (DS)  component,  Q  to  the  extended  source  (ES) 
component  (in  this  case  colored  noise),  and  IN  (an  N  by  N  identity  matrix)  to  the  component  due 

to  uncorrelated  noise.  Normalize  P  and  Q  so  that  each  trace  equals  N.  The  three  variances, 
<&.  <4.  and  Oq,  represent  the  power  at  a  single  sensor  due  to  discrete  sources,  extended 

sources,  and  uncarrelated  noise,  respectively. 

Define  the  discrete  source  component  P  corresponding  to  q  £  N  discrete  angular  sources 
at  Oj,  i  =  1,  2, ....  q,  as^ 

P  =  DCDh  ,  (4) 

where  the  N  by  q  steering  matrix  D  =  [dj  d2  ...  dq  ]  contains  die  i-th  column 

d f  =  [l  e'*,d  e'ik,<24,> ...  e‘*,(N'1>d] ,  (5) 

and  die  wavenumber  of  the  i-th  source  at  is 

ki  =  ycosft) .  (6) 
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The  factor  d  is  the  nominal  spacing  between  sensors,  given  as  a  fraction  of  the  wavelength  X. 
Write  the  above  expression  as  a  function  of  angular  frequency  equal  to  the  cosine  of  the  bearing 
angle  of  interest  6,  with  0  £  0  £  n.  Assume  a  diagonal  source  coherence  matrix  C  so  that  no 
correlation  exists  between  sources.  A  thin  spike  or  delta  function  illustrates  the  spatial  density  of  a 
discrete  source  at  bearing  angle  6. 

EXAMPLE  1:  EXTENDED  SOURCE  ENERGY  SPECTRUM 

The  angular  energy  spectrum  of  an  extended  source  at  bearing  angle  6  affects  a  spatial 
bandwidth  around  the  angle  6. 

Compute  the  spatial  density  of  an  arbitrary  extended  source  modeled  as  a  first  order 
autoregressive  (AR)  process  with  coefficient  a,  =  -0.7  steered  to  45  degrees.  Figure  1  depicts  the 
spatial  density. 


Figure  1.  Extended  Source  Spectrum 

Model  the  P  extended  sources  in  component  Q  as  the  sum  of  L-th  order  AR  spatial 
processes  Q;  utilizing  the  Gohberg  Formula^  with  complex  coefficients  a;  j,  scaled  by  a  factor 

otp  i  =  1,  2, ...,  P,  as  follows: 
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(7) 


(8) 

(9) 


A  useful  description  of  a  complex  AR  parameter  is  the  polar  form,  consisting  of  magnitude 
and  steering  angle.  Given  the  magnitude  of  an  AR  coefficient,  a,,  and  a  steering  angle,  6, 
compute  the  equivalent  complex  coefficient  a,  as 

a,  =  a,ej2*dco*(9) .  (11) 

Later  in  this  report  when  the  standardized  test  case  (STC)9  is  used  to  test  the  prewhitening 
algorithm,  the  expression  in  equation  (1 1)  will  be  useful. 

EXAMPLE  2:  COMPLEX  AR  COEFFICIENT 


For  an  array  with  quarter-wavelength  sensor  spacing  (d  =  0.25),  model  the  extended 
source  from  example  1  by  computing  the  equivalent  complex  AR  coefficient  as  follows: 

a,  =  -0.7,  (12) 


5 


8  =  45  degrees  , 


(13) 


and 

a,  =  -0.7ej2*(0‘25)o“(45) 

=  -0.3108  -  0.62721.  1 

While  figure  1  has  illustrated  the  spatial  energy  density  for  this  process,  figure  2  shows  a  pole  plot 
of  the  coefficient  a,. 
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Figure  2.  Pole  Pkr  ‘or  Example  2 

EXAMPLE  3:  MUSIC  IN  A  COLORED  BACKGROUND  NOISE 

Now,  consider  an  example  that  illustrates  the  degradation  in  performance  when  the  MUSIC 
direction  finder  is  used  in  a  nonwhite  background  spectrum  environment 

Given  an  eight  sensor  array  with  quarter-wavelength  sensor  spacing  and  a  0-dB  extended 
source  from  example  2,  inject  two  discrete  sources  into  the  observed  noise  field: 


6 


Discrete  source  1:  30  degrees,  SNR  =  -3  dB 
Discrete  source  2:  70  degrees,  SNR  =  -19  dB. 

Compare  the  MUSIC  spectra  for  the  prewhitened  and  non-prewhitened  cases.  For  the  MUSIC 
computation,  assume  the  presence  of  q  =  3  sources. 

In  figure  3,  the  dotted  vertical  lines  represent  the  true  direction  of  arrival  (DO A)  angles. 
The  dashed  line  corresponds  to  the  pre whitened  MUSIC  spectra.  The  non-prewhitened  case 
exhibits  a  strong  bias  due  to  the  extended  source  at  45  degrees.  The  stronger  extended  source 
obscures  the  low  SNR  source  at  70  degrees. 


Figure  3.  MUSIC  Spectra  for  Example  3 
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METHOD  OF  PREWHITENING 


Define  x  as  an  N-dimensional  vector  whose  covariance  matrix  is  given  by 
R  =  E[xxH] 

=  R-<#n  (15) 

where  x  =  [x,  x2  ...  xN]T.  Note  that  R  is  the  covariance  matrix  in  the  absence  of 
uncorrelated  noise  and  can  be  estimated  by  using  the  smallest  eigenvalue  of  R  as  an  estimate  of  a„. 

If  Les  is  the  lower  triangular  matrix  Cholesky  factor  given  by 

LesLes  =  OesQ  ,  (16) 

defme  a  prewhitened  snapshot  vector  as 

xpw  =  L^x  .  (17) 

The  covariance  matrix  of  the  prewhitened  vector  then  contains  a  purely  white  spatial  noise 
density. 


Rpw  =  eJx^x,^,] 

=  L'^RL” 

=  ^  ds  Ees  PEes  +  OgLjjQLgj 

=  a^PLS  +  <414 


ES 


'ES 


—  rr2  T'1  PT'H  +  I 

—  ^ds^es'^es  t  *n  » 


where 


(18) 


(19) 
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The  algorithm  for  prewhitening  the  measured  matrix  R  minimizes  a  likelihood  functional  * 1  to 
obtain  an  estimate  of  Q. 

THE  LIKELIHOOD  FUNCTIONAL 


Following  the  development  in  reference  (11),  assume  we  are  given  K  array  data  snapshots 
in  a  sequence, 

{x,,  x2,  xK}  ,  (20) 

with  covariance  R.  Form  the  estimate1 2  of  R  as 


i 

-  Yxkx?  . 
K  ^  k  k 


k*0 


(21) 


A 

The  likelihood  functional  relates  the  estimate  R  to  an  assumed  covariance  structure  R.  To  aid  in 
the  derivation,  first  assume  that  the  true  covariance 


R  »  P  +  XIN 


(22) 


corresponds  to  the  white  noise  case,  and  omit  the  scale  factors  for  brevity.  The  factor  X  equals  a 
positive  scalar  parameter.  The  matrix  P  equals  an  unknown  covariance  matrix  of  rank  q. 
Although  methods1 3- 15  exjst  for  the  computation  of  q,  the  number  of  sources  present  assume  that 
q  represents  a  known  quantity.  The  function 

f({x„  x2,  ...,  xK}  /  X,  P)  =  ffKW|  k  e-*«~(RK  )  (23) 

represents  the  conditional  density  for  the  complex  Gaussian  sequence  xjf  i  =  1,  2,  ...,  K,  with 
covariance  R  conditioned  on  (X,  P)  and  1 1  being  the  determinant  operator.  The  function  trace  (A) 
sums  the  diagonal  elements  of  matrix  A.  The  problem  now  involves  seeking  the  values  for  (X,  P) 
that  maximize  this  likelihood.  To  express  the  log-likelihood  as  a  functional  in  terms  of  the 

A  A  A  A 

eigenvalues  of  R ,  X,  >  X2  >  ...  >  XN,  define 
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Wq>  =  n  x  (M) 

and 

,  (  N  \ 

Oq)  =  rp-  I  \  (25) 

N-q  J 

as  the  geometric  and  arithmetic  means,  respectively,  of  the  N-q  lowest  eigenvalues  of  R .  As 
shown  in  reference  (1 1),  we  can  then  write 

max  ln(f({x,,  x2,  ....  xK}  /  X,  P)) 

=  K  f-N  In  (Jt)  -  N  -  In  (|r|)  +  £ln  (*i)  -  (N-q)  In  (Wq))!  .  (26) 

<  i»q+l  j 

Maximizing  the  likelihood  is  equivalent  to  maximizing  die  functional: 

L,(Q)  =  (N  -q)(ta  (XGM(q))  -  In  (XAM<q)))  .  (27) 

For  the  spatially  nonwhite  case,  use  the  form 

R  =  P  +  XQ  (28) 

for  the  assumed  true  covariance  matrix.  As  in  the  preceding  section,  decompose  Q  into  Cholesky 
triangular  factors, 

Q  =  LLh  ,  (29) 

and  form  the  prewhitened  covariance  matrix 

Rpw  =  L'RL"  .  (30) 
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It  can  then  be  shown  (see  reference  (11))  that  the  likelihood  functional  conditional  to  the  noise 
covariance  Q  with  q  sources  is 

MQ)  =  (N-q)(ln  (l  WGM  <q))  -  In  (X,WAM(q  ))) ,  (31) 

where  A.WGM(q)  and  XWAM(q)  correspond  to  the  geometric  and  arithmetic  means,  respectively,  of 
the  eigenvalues  of  the  prewhitened  covariance  matrix  .  The  geometric  and  arithmetic  means 

realize  equality  only  for  the  case  of  equality  between  the  N-q  lowest  eigenvalues.  The  equality  of 
the  N-q  smallest  eigenvalues  represents  a  flat,  spatially  white  noise  spectrum,  and  in  this  case, 
Lq(Q)  equals  zero. 

MAXIMIZATION  OF  THE  LIKELIHOOD  FUNCTIONAL  -  LE  CADRE'S 
TECHNIQUE 

The  prewhitening  problem  next  becomes  one  of  maximizing  Lq(Q)  relative  to  the 
coefficients  a,,  a2,  ....  aL  parameterizing  the  noise  covariance  matrix  Q.  Accomplish  the 
maximization  by  an  iterative  gradient  algorithm1 1  at  step  k  with  parameter  vector 

A;  =  (o2k,  a*,  a2k,  ...,  ak)  (32) 

and  general  form 

Ak+i  =  Ak  -  p0Gk ,  (33) 

where  p0  defines  a  step  size,  in  this  case  constant,  and  Gk  defines  the  gradient  vector. 
Computation  of  the  i-th  element  of  Gk  consists  of  calculating 

c,(i)  =  (34) 

for  i  =  1,  2, ...,  L.  Figure  4  expresses  a  functional  relationship  between  the  eigenvalues  of  the 
prewhitened  covariance  matrix,  the  likelihood  function,  and  the  AR  coefficients.  The  algorithm 
seeks  to  find  the  maximum  of  the  likelihood  functional  surface  parameterized  by  the  AR 
coefficients. 
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Figure  4.  Relation  of  Likelihood  Functional  to  ARMA  parameters 

EXAMPLE  4:  LIKELIHOOD  FUNCTIONAL  SURFACE  CONTOUR 

Given  the  acoustic  field  from  example  3,  plot  the  maximization  surface  contour  of  Lq(Q) 
in  figure  5  as  a  function  of  the  real  and  imaginary  parts  of  the  coefficient  a, . 


Real 


Figure  5.  Maximization  Surface  Contour  Plot  of  Lq(Q) 
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Clearly,  the  maximum  value  of  Lq(Q)  occurs  at  the  exact  value  of  the  coefficient  computed  in 
example  2. 


Computation  of  the  gradient  vector  reduces  to 

G„(i)  -  ^  (35) 


for  i  »  1,  2, ....  L  and  j  =  1,  2, ....  N. 

LEMMAS:  EIGENVALUES  OF  PREWHITENED  COVARIANCE  MATRIX 


The  following  lemmas  represent  intermediate  steps  in  this  computation. 

Lemma  1:  The  eigenvalues  of  Rw  k  equal  the  eigenvalues  of  QkR. 

Proof:  To  compute  the  eigenvalues  X,,  X2,  ....  XN  of  the  N  by  N  matrix  A,  solve 

|A  -  UN|  -  0. 

With  this  fact  and  defining 
LkLk  =  Qk  , 

then 

|Qk'R  -  uN| 

=  |LkLHkR  -  XIN| 

=  |LHkRLk  -  XINj 

=  I^PW.  k  ■  ^n|  • 


(36) 


(37) 


(38) 


Lemma  2:  If  R  =  TTH,  the  eigenvalu  *■  R  V  k  equal  the  eigenvalues  of  THQkT. 
Proof: 
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|q,'r  -  u,| 

=  (q,’tth  -  u„| 
=  frMQ;'T  -  xiH| . 


(39) 


Therefore  the  computation  of  dXj  /  daf  reduces  to  the  computation  of  dp*  /  da* ,  the  partial 
derivatives  of  the  eigenvalues  of  the  Hermitian  matrix  THQ'kT .  Since  the  AR  coefficients  form  an 
explicit  function  for  Qk ,  this  reduction  represents  an  important  result^  1  Additionally,  use  the 
following  classical  result  for  an  N  by  N  Hermitian  matrix  A  (A  =  AH)  with  simple  eigenvalues 
Xj  and  associated  eigenvectors 


•  3  A 

v'a7Av. 


(40) 


Equation  (40)  derives  from  the  eigenstructure  definition 


A 


j-i 


(41) 


ALGORITHM  FOR  MAXIMIZATION  OF  THE  LIKELIHOOD  FUNCTIONAL 
BY  COMPUTATION  OF  THE  GRADIENT  VECTOR  Gk 


Continuing  to  parallel  the  development  in  reference  (1 1),  we  employ  the  above  definitions 
and  lemmas  to  write  the  algorithm  as  follows. 

1 .  Compute  partial  derivatives  of  Qk  relative  to  the  complex  AR  parameters  a, : 


for  i  as  0,  1, ....  L. 

(See  appendix  A  for  derivation  of  this  expression  and  the  definition  of  Z'.) 
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S.  Select  step  size  to  ensure  adequate  convergence  rate. 

Although  methods^  1  exjst  for  updating  p0,  assume  a  fixed  step  size  arrived  at  by  trial-and-error. 
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PERFORMANCE  METRICS 


Four  performance  metrics  will  be  used  to  evaluate  algorithm  performance: 

1.  Compute  a  pole  plot  of  the  AR  coefficients  by  generating  one  at  each  iteration  of  gradient 
algorithm  to  view  convergence  trajectory. 

2.  Compute  the  cosine  of  the  angle  between  the  estimated  covariance  and  the  true  covariance 
Qtrue- 


COs(Qe3T*  QtRUe)  = 


trace(QESTQTRUE) 

V  ^^(QestQestJ^^'^QtrueQtrue  ) 


(48) 


A  cosine  close  to  1  implies  collinearity. 


3.  Compute  the  AR  spatial  density  for  estimated  coefficients: 


S(z) 


1 

A(z)A*(z)  ’ 


where  "  * "  corresponds  to  the  complex  conjugate  and 
A(z)  =  a0  +  a,z  +  a2z2  +  ...  +  aLzL 
z  =  c*i2,rtC0,(#) 

4.  Compute  the  MUSIC  26  bearing  response  for  prewhitened  case: 


MUSIC 


(8)  = 


1 


(  N  ’N 

eH(0) 


k(0) 


(49) 


(50) 

(51) 


(52) 


where 


eT(0)  =  fl  ej2«*«(«)  cj4»*«(6) 


ej2(N-l)itdra»(9) 


(53) 
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and  u  j  represents  an  eigenvector  of  the  prewhitened  covariance  matrix  Up* . 


THE  STC  SIMULATIONAL  DATA 


To  evaluate  the  effectiveness  of  the  proposed  method,  use  the  STC,^  a  realistic  model  with 
a  challenging  range  of  SNR  and  difficult  source  locations.  The  STC  employs  quarter-wavelength 
sensor  spacing.  Figure  6  depicts  the  STC  source  locations  graphically  and  table  1  enumerates  the 
relative  source  strengths. 


Cos  (Bearing) 
Figure  6.  The  STC 


Table  1.  The  STC 


COS(6) 

AR 

PARAMETER 

SNR 

(dB) 

DISCRETE 

SOURCES 

1 

105 

-026 

- 

-20 

2 

90 

0 

- 

-3 

3 

85 

0.09 

- 

-3 

4 

70 

0.34 

- 

0 

5 

45 

0.71 

- 

-14 

EXTENDED 

SOURCES 

1 

0 

1 

-0.7j 

0 

2 

135 

-0.71 

-0.4263+0.8602i 

-27 
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RESULTS  OF  SIMULATION 


Assume  an  L  =  4  AR  model  and  initialize  the  algorithm  with  parameter  vector 


At  =  (1,  0,  0,  0,  0) .  (54) 

Run  the  algorithm  for  25  iterations  with  step  size  p0  =  >0.025.  Set  q  equal  to  the  number  of  signal 
eigenvalues,  7,  for  the  prewhitener  and  for  MUSIC.  Assume  a  known  noise  floor  ojj. 

Figures  7  through  10  illustrate  the  results  of  the  Le  Cadre  algorithm  chi  the  STC.  The  great 
disparity  in  SNR  for  the  two  extended  sources  (ES)  causes  the  algorithm  to  prewhiten  only  the 
noise  associated  with  ES  1  at  0  degrees.  Although  a  fourth-order  AR  process  attempts  to  model 
the  sum  of  two  first  order  processes,  in  this  case,  a  first-order  model  suffices.  The  MUSIC 
spectrum  shows  that  although  the  prewhitening  of  ES  1  allows  detection  of  the  discrete  source  at 
45  degrees,  it  also  causes  the  promotion  of  a  weak  second  extended  source,  ES  2.  The  number  of 
assumed  sources  plays  an  important  role  in  the  detection  of  ES  2;  indeed  if  q  s  5,  the  peak  at  135 
degrees  disappears.  Attempts  to  prewhiten  the  small  source  at  135  degrees  with  a  second  pass  of 
the  algorithm  fail  because  of  the  nature  of  the  remaining  eigenvalues  of  the  prewhitened  R. 
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Cosine  Metric 


1 


0.5 


0 


-0.5 


-U - . - - — - - • - 1 

-1  -0.5  0  0.5  1 

Real 


Figure  7.  Pole  Plot  for  STC 


Figure  8.  Cosine  Metric  for  STC 
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cos(bearing) 


Figure  9.  Noise  Spatial  Density  for  STC  Figure  10.  MUSIC  Response  for  STC 

Figure  1 1  depicts  the  eigenvalues  associated  with  the  STC.  The  solid  line  corresponds  to 
the  non-prewhitened  case  and  the  dotted  line  corresponds  to  the  case  following  the  25th  iteration  of 
the  algorithm.  The  dashed  line  represents  the  eigenvalues  omitting  the  two  extended  sources 
entirely.  To  visualize  the  failure  of  the  attempted  second  pass  of  the  algorithm,  recall  that  the 
likelihood  functional  measures  the  proximity  of  the  N-q  (24  -  7  =  17)  smallest  eigenvalues.  In 
this  case,  a  comparison  of  these  eigenvalues  illustrates  the  algorithm  difficulty.  The  dashed  line 
lies  so  close  to  the  dotted  line  that  the  gradients  involved  in  the  computation  of  the  update  become 
small  relative  to  the  step  size.  A  flat  eigenvalue  spectrum  corresponds  to  a  value  of  zero  for  the 
likelihood  functional. 
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Figure  11.  Eigenvalues  of  the  STC 

Table  2  shows  the  values  of  the  likelihood  functional  for  the  various  covariance  structures. 
The  magnitude  of  the  likelihood  functional  directly  relates  to  the  size  of  the  gradient  and 
consequently  to  the  convergence  rate.  While  the  first  pass  converged  in  about  10  iterations, 
attempts  at  modifying  the  step  size  for  the  second  pass  resulted  in  extremely  slow  convergence 
and,  ultimately,  algorithm  failure. 

Table  2.  Values  of  Likelihood  Functional  for  STC 


Lq(Q)  =  (N-q)(ln(XGM(q))  -  In^Cq))) 

R  =  als  P  +  olsQ  +  <JoIN 
(solid) 

-0.5501 

Rw  =  LRLh 
(dotted) 

-0.0009442 

(dashed) 

0 

To  visualize  the  gradients  involved,  figure  12  depicts  the  likelihood  functional  of  the  STC 
as  a  function  of  the  real  and  imaginary  components  of  the  primary  AR  coefficient.  In  this  case,  the 
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primary  pole  equals  -0.0044  -  0.6910i,  the  maximum  of  the  surface.  Figure  13  represents  a 
contour  of  the  surface. 


Figure  12.  Likelihood  Functional  Surface  for  STC 


Real 


Figure  13.  Likelihood  Functional  Contour  for  STC 

A  second  interpretation  of  the  failure  of  the  second  pass  of  the  algorithm  relies  on  the 
magnitude  of  the  AR  coefficient  for  the  second  extended  source.  The  magnitude  of  the  coefficient 
equals  0.96,  close  to  the  unit  circle.  As  the  AR  coefficient  moves  closer  to  the  unit  circle,  the 
Gohberg  formula  produces  an  extended  source  covariance  with  extreme  ill-conditioning  due  to  a 
very  strong  principal  eigenvalue.  Essentially  this  type  of  covariance  models  a  discrete  source.  The 
source  contributes  only  to  the  principal  eigenvalues  of  the  covariance  R.  In  this  manner,  the  pre¬ 
whitening  eigenvalues  that  appear  after  the  principal  eigenvalues  remain  unaffected,  contributing  to 
a  flat  white  eigenvalue  spectrum. 
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CONCLUSIONS 


The  STC  emphasizes  certain  deficiencies  of  the  Le  Cadre1  1  prewhitening  technique.*  The 
presence  of  multiple  extended  sources  that  comprise  the  spatially  colored  noise  results  in  the 
greatest  difficulty.  The  great  disparity  in  SNR  between  the  two  extended  sources  causes  the 
algorithm  to  concentrate  almost  solely  on  the  stronger  source.  After  prewhitening  the  first  strong 
source,  a  promotion  of  the  weaker  source  occurs  using  MUSIC,  although  the  noise  eigenvalues 
represent  nearly  a  white  noise  spectrum.  A  second  pass  of  the  prewhitener  fails  due  to  the  nature 
of  the  prewhitened  eigenvalues.  Although  the  results  are  not  presented  here,  inaccurate  knowledge 
of  the  noise  floor  diminishes  algorithm  performance  substantially. 

The  estimate  of  the  number  of  sources  present  in  the  acoustic  field  also  greatly  affects 
algorithm  performance.  Since  many  problems  in  sensor  array  processing  rely  on  this  estimate, 
research  into  its  computation  actively  continues. 

All  optimization  problems  based  on  gradient  algorithms  suffer  from  inaccurate  initialization 
and  computation  of  suitable  step  size.  Indeed,  Le  Cadre1 1  presents  a  method  for  step  size 
updating  and  provides  the  algorithm  starting  point  Other  suitable  techniques  for  optimization  exist 
for  this  type  of  problem,  notably  the  Broyden-Fletcher-Goldfarb-Shanno  routines.^? 

The  parametric  assumptions  on  the  noise  model  restrict  the  utility  of  this  approach  to 
nonuniform  arrays.  The  AR  formulation  fails  when  inaccurate  sensor  spacing  or  a  curvature  exists 
in  the  line  array.  Future  work  will  attempt  to  address  these  inadequacies.  Recent  related  research 
includes  the  work  of  Kay  and  Nagesha.^8, 29 


*  The  MATLAB  source  code  in  this  report  will  be  provided  to  any  interested  party  by  electronic  media. 
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APPENDIX  A:  COMPUTATION  OF 

da; 


(Omit  all  unnecessary  indexing  for  simplicity.) 
Define  the  N  by  N  matrices  Z1  such  that 


J 1  ifj-k  =  i  1  £  j  £  N 
{o  else  1  <  k  <  N 


Write  the  Toeplitz  matrices  A,  and  A3  as 


i  «  1 


A,  =  £a,Z'  =  a0IN  +  . 

i  -  0 

A,  =  £a,Z" ' . 


To  simplify  the  computation,  formulate  A,  and  A3  in  vector  form: 


A|  —  +  [ai  a2  •”  aL 


Z‘ 

Z2 

ZL 


=  a0IN  +  aTZ, , 


A3  -  [a,  a2 ...  aL] 


zN-' 

zN-2 

r»N-L 


=  a  Z2  , 


(A-l) 


(A-2) 


(A-3) 


(A-4) 


where  emphasizes  a  vector  variable.  Express  the  complex  coefficients  a,  as  the  sum  of  real 
and  imaginary  parts  x{  and  y{  so  that 


A-l 


a<>  *  *o  +  jy0 . 

a  =  x  +  jy  , 

where  j  =  V-I.  With  the  Gohberg  formula  formulate  the  matrix  Q  as 
Q '  =  J-  (A, A?  -  A, A,") 

°ES 

i  f{(x0  +  jy0)iN  +  (*T  +  jyT)z,}{(x0  -  jy0)iN  +  %(*  -  jy)} 

-  {(xT  +  jyT)Z2}{zJ(x  -  jy)} 

First  compute  the  partial  derivatives  with  respect  to  the  zero-th  coefficient: 
dQ-'  _  l  !n{(xo  -  Jyo)IN  +  Zf(x  -  jy)} 

3x°  +  {(x0  +  jyo)iN  +  (*T  +  jyT)z,  }iN 

=  4-[2x0In  +  zj(x  -  jy)  +  (xT  +  jyT)Z,] , 

=  ~-[A,  +  A”] , 


i  jiN{(x0  -  jy0)iN  +  -  jy)} 

+  {(x0  +  jy0)lN  +  (xT  +  jyT)Z,}(-jIN) 

j_  (x0  -  jy0)*N  +  %(*  -  jy)  -  (xo  +  jyo)iN 
°es  -  (xT  +  jyT)Z, 

-  -i-W  -  A-J 


(A-5) 


(A-6) 


(A-7) 


(A-8) 


Define  the  complex  gradient^  as 


aol  .  .aol 

v.  ^  J  -v. 


-  ^  ■ 
®ES 


(A-9) 


A-2 


Compute  the  partial  derivatives  with  respect  to  the  first  through  the  Lth  coefficients: 

z,{(*.  -  jy«)I«  +  Z,T(*  -  jy)} 

+  {(*o  +  jy»)iN  +  (*T  +  jyT)z,}z,  . 

-  z,(z;(x  -  jy))  -  (iiT  +  jyT)Z,(Z,) 

=  4-[Z,Ar  +  A,Z,  -  Z,A,"  ■  AjZj]  , 

C^es 

=  4-[jZ,Ar  -jA,Z,  -  jZ.A?  +  jA)Zj] , 
oy  Oes 

=  4-fz.Ar  -  A,Z,  -  Z,a;  +  A3Zj]  , 

ctes 

dQ'  3Q1  dQ' 
da  ~  diT  JdyT  ’ 

=  ^-[a.Z)  -  A,ZJ] . 

S 


dQ  '  _  _1_ 
dx  Ob 


(A- 10) 


(A-ll) 


(A-12) 


The  expressions  in  equations  (A- 10)  -  (A- 13)  also  correctly  compute  the  partial  of  the  zero-th 
coefficient  since 

Z°  =  IN  (A-13) 

and 

ZN  =  0N  (A- 14) 
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APPENDIX  B:  SOURCE  CODE  FOR  PREWHITENING  ALGORITHM 


%  ******************************************************* 
%  This  MATLAB  program  implements  Le  Cadre's  algorithm  for 
%  estimating  the  noise  covariance  matrix  B  for  the  STC. 

%  Author:  Alain  C.  Barthelemy 
%  Date:  23  June  1993 

%  ******************************************************* 


clear; 

clg; 

clc; 

disp(' '); 
hold  off; 
axis  ( 'normal' )  ; 


%  PARAMETERS  FOR  STANDARDIZED  TEST  CASE  (STC) 


0 

th 

vs 

vamoiz 

varsig 

ambnoiz 

sp 

N 

thn 


sqrt  (-1) ;  %  Cooplex  variable 

[45  70  85  90  105] ;  %  Discrete  source  bearings 

[0.02  0.515  0.23  0.23  0.005];  %  Source  variances 


lj 
2; 

1; 

0.25; 

24; 

[0  135] 


%  Fewer  of  noise 
%  Fewer  of  signal 
%  Fewer  in  ambient  noise 
%  Nominal  sensor  spacing 
%  Total  Number  of  sensors 
%  DQA  angles  for  extended  sources 


%  Complex  AR  coefficients 

a  =  [-0.7  -0.96]  :*e>?5(j*2*pi*sp*cos(t±n*pi/180) )  ; 


alpha  =  [0.998  0.002]; 


%  Weighting  for  2  extended  sources 


%  PARAMETERS  FOR  PRE-WHITENING  ALGORITHM 
p  =  4;  %  Estimated  order  of  AR  process 

q  =  7;  %  Estimate  of  the  number  of  sources 

rho  =  -0.025;  %  Step  size  for  gradient  algorithm 

num_k  =  25;  %  Maximum  number  of  Le  Cadre  iterations 


%  PARAMETERS  FOR  DISPLAY 


%  [#  circle  divs,  #  polar  lines,  #  circles] 


B-l 


grjp  =  [100  12  4]; 

theta  =  linspace(0,  2*pi,  gr_p(l)); 

phi  =  linspace(0,  2*pi-(  (2*pi) /gr_p(2) ) ,  gr_p(2)); 

this  =  cnes(gr_p(l)  ,1)  *linspace(l/gr_p(3)  ,l,gr_p(3) )  ; 


%  COMPUTE  COVARIANCE  STRUCTURE  FOR  STC 


B01  =  B_F(vamoiz,  a(l),  N)  ; 

B01  =  (N/trace(B01))*B01; 

B02  =  B_F(vamoiz,  a  (2),  N)  ; 

B02  =  (N/ trace (B02 ))*B02; 

BO  =  alpha  (1)  *B01  +  alpha(2)  *B02; 
%  5  discrete  sources 


%  First  extended  source 
%  Normalize 

%  Second  extended  source 
%  Normalize 
%  Scaled  sum 


DP  =  exp(sqrt(-l)*2*pi*sp*[0:N-l]  ’*(cos(th*pi/180) ) )  ; 

P  =  DP*diag(vs)*DP';  %  Scaled 
P  =  (N/trace(P)  )*P;  %  Normalize 
R  =  varsig*P  +  vamoiz*B0  +  airfcnoiz*eye(N) ;  %  Covariance 


%  REMOVE  DIAGONAL  FROM  R,  ASSUME  KNOWLEDGE  OF  NOISE  FLOOR 


R  =  R  -  amtxioiz*Qre  (N)  ; 

%  R  NOW  CONTAINS  SMALL  IMAGINARY  COMPONENTS  ON  ITS 
%  DIAGONAL  MAKE  R  POSITIVE  DEFINITE  FOR  CHOLESKY 
%  COMPUTATION,  THIS  DOES  NOT  ALTER  R  IN  ANY  SIGNIFICANT 
%  WAY 


R  =  (R+R')/2; 

%  COMPUTE  A  YULE-WALKER  APPROXIMATION  OF  THE  POLE 
%  LOCATIONS,  OF  THE  SUM  OF  THE  TWO  EXTENDED  SOURCES 

poles  =  B0(l:p,  l:p)  \  ( -BO  ( 2 :  (pfl) ,  1)); 

B_poles  =  B_FI  (vamoiz,  poles,  N) ;  %  Use  in  cosine  metric 

%  PRE-COMPUTE  Z  MATRICES 

Z1  =  zeros  (  (ph-1)*N,  N) ;  %  Initialize 

Z2  =  Zl; 
for  i  =  0:p, 

Zl(  (1:N)  +i*N,  :)  =  Z_FDNC(i,  N)  .  ' ; 

Z2  ( (1  :N)  +i*N,  : )  =  Z_FTJNC (N-i,  N) . '  ; 
end; 
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%  LB  CADRE'S  ALGORITHM  FOR  ESTIMATING  THE  COVARIANCE 
%  MATRIX  B 


A  =  [  1 ;  zeros  (p,  1)];  %  Initialize  parameter  vector 

T  =  chol  (R) ' ;  %  Caqpute  Cholesky  factor  of  R 

metric  =  zeros  (num_k,  1);  %  Initialize  cosine  metric 

G  =  zeros  (p+1,  1);  %  initialize  gradient 

%  Initialize  matrix  of  estimated  parameters 
A_plot  =  zeros (nuiLK,  (p+1)); 


for  k  =  l:nurn_k. 


%  Iteration  index 


A  =  A  -  rho*G;  %  Update  parameter  vector 

A_plot  (k,  : )  =  A. ' ;%  Fill  in  matrix  of  estimated  parameters 
%  Form  estimated  covariance  matrix  B 

invJB  =  B_FT  (A(l) ,  A(2 :  (p+1) ) ,  N)  ; 


%  COMPUTE  PERFORMANCE  METRIC 

metric (k)  =  ... 

real  (trace (inv_B*B_poles)  /sqrt  (trace (B_poles*B_poles)  *  . . . 
trace  ( inv_B*inv_3) ) )  ; 

%  CALCULATION  OF  GRADIENT  VECTOR  G  REQUIRES  KNOWLEDGE  OF 
%  THE  EIGENSTRUCTURE  OF  T'*INV(B)*T 

[U,  LAM]  =  eig(T'*inv_B*T) ;  %  C3crapute  eigenstructure 

%  Vector  of  small  eigenvalues 

LAM  =  di.ag(LAM(q+l:N,  q+l:N)); 

ar  =  mean  (LAM) ;%  Arithmetic  mean  of  small  eigenvalues 

%  UPDATE  OF  GRADIENT  VECTOR  G 


for  i  =  0:p,  %  Index  over  AR  coefficients 

%  PARTIAL  DERIVATIVES  OF  INV(B)  RELATIVE  TO  A  SUB  I 

deli  =  (2/A(l))*(Al_F(A(2:  (p+1)),  N)*Z1(  (l:N)+i*N,  :)  ... 

-  A3_F(A(2:  (p+1)),  N)  *Z2  ( (1  :N)  +i*N,  :))? 

%  DERIVATIVE  MATRICES 
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deliprime  =  T'*deli*T; 

%  PARTIAL  DERIVATIVES  OF  THE  EIGENVALUES  OF  PRE-WHITENED  R 

partial  =  diag  (U  ( : ,  q+1  :N)  '  *deliprime*U  ( : ,  q+1  :N) )  ; 

%  UPDATE  OF  GRADIENT  VECTOR  G 

G(i+1)  =  sum  (partial. /LAM)  -  sum  (partial  )/ar; 
end; 
end; 

%  SET  UP  DISPLAY  GRAPH  AND  POLE  PLOT 

axis  ( '  square ' ) ;  %  Set  aspect  ratio  to  square 

axis([-l  1-11]);  %  Set  coordinate  axes 

%  PLOT  POLE  PLOT 

plot  (real  (poles) ,  imag  (poles ) ,  'clo',  [zeros(grjp(2) ,  1)  ... 
cos(phi)']',  [zeros(gr_p(2) ,  1)  sin(Fhi)']',  'cl:',  ... 
this.* (ones (gr_p(3) ,l)*cos (theta) )’ ,  this.*  ... 

(ones (gr_p (3 ) ,l)*sin( theta) )' ,  'cl:',  ... 

real (A_plot ( : ,  2:(pfl))),  imag (A_plot ( : ,  2:(p+l))),  'clx'); 

xlabel  ( 'Real ' ) ;  %  Label  x  axis 

ylabel  ( '  Imaginary ' ) ;  %  Label  y  axis 

pause;  %  Pause  in  execution 

%  PLOT  COSINE  METRIC 

axis  ( 'normal ' ) ;  %  Set  aspect  ratio  to  default 

%  Size  coordinate  axes 

axis  ( [1  nurn_k  min  (metric)  max  (metric)  ] )  ; 

%  PLOT  METRIC 

plot ( 1 :num_k,  metric,  ' cl- ' ) ; 
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grid; 

xlabel  ( '  Iteration  number ' )  ; 
ylabel  ( 'Cosine  metric' ) ; 


%  Put  grid  onto  graph 
%  Label  x  axis 
%  Label  y  axis 

pause;  %  Pause  in  execution 

%  COMPUTE  AND  PLOT  NOISE  SPATIAL  DENSITIES 

M  =  301;  %  Number  of  bearing  h-inm 

%  Angles  from  0  to  180  degrees 

degree_angles  =  linspace(0,  pi,  M)  ; 

bearing  angles  =  cos  (degree_angles) ;  %  Cosine  of  above 

psd_plot  =  zeros  (M,  num_k) ;  %  Initialize  far  plotting 

for  k  =  l:num_k,  %  iteration  index 

%  COMPUTE  NOISE  SPATIAL  DENSITY 

psdjplot  ( ; ,  k)  =  ARMA([1  A_plot(k,  2:(p+l))],  1,  ... 
degree_angles,  sp) . ' ; 

%  ADD  ARBITRARY  FACTOR  FOR  DISPLAY  PURPOSES 

psd_plot(:,  k)  =  psdjplot  ( : ,  k)  +  (k*3); 
and; 

axis([l  2  3  4]);  axis;  %  Set  axes  scaling  to  automatic 

%  PLOT  NOISE  SPATIAL  DENSITY 

plot(bearingL.angles,  psd_plot,  'cl-'); 

xlabel  ( '  cos  (bearing) ' ) ;  %  Label  x  axis 

ylabel  ( '  Iteration  ->' ) ;  %  Label  y  axis 

pause;  %  Pause  in  execution 

%  COMPUTE  AND  PLOT  MUSIC  RESPONSES 

e  =  le+16;  %  Enhancement  factor  for  signal  subspace 

DS  =  e>¥>(-sqrt(-l)*2*pi*sp*(0:N-l) . ' . . . 

*cos  ( 0 :  (pi/(M-l) )  :pi) ) ;  %  Steering  matrix 
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num_src  =  q;  %  Assumed  number  o£  sources  £or  MUSIC 

br_plot  =  zeros  (M,  num_k) ;  %  Initialize  bearing  response  plot 

for  k  =  l:num_k,  %  Iteration  index 

%  FORM  ESTIMATED  NOISE  COVARIANCE  MATRIX  B 

B_EST  =  B_F  (Ajplot  (k,  1) ,  A_plot(k,  2:  (p*l) ) .  ' ,  N)  ; 

B_EST  =  (N/trace  (B_EST) )  *B_EST;  %  Normalize 

%  Make  positive  definite  for  Cholesky 
B_EST  =  (B_EST+B_EST ' ) /2; 

L  =  chol  (B_EST) ' ;  %  Caspute  Cholesky  factor 

%  FORM  PRE-WHITENED  ESTIMATE  OF  R 

RW  =  inv(L)*R*inv(L' )  ; 

%  ADD  NOISE  FLOOR  BACK  IN 

RW  =  RW  +  ambnoiz*e/e  (N)  ; 

[U,  LAM,  AA]  =  svd(RW) ;  %  Canpute  eigenstructure  of  RW 

LAM  =  diag(LAM) ;  %  Vector  of  eigenvalues 

%  COMPUTE  MUSIC  BEARING  RESPONSE 

br^plot ( : ,  k)  =  ... 

real  ( ( [sum(DS. *  ( ( (eye (N)  -U( : ,  (1  :num_src) )  * . . . 
diag(e*(LAM(l:num_src)-l)  ./(l+e*(LAM(l:nuin_src)-l) ) ) . . . 

*U( : ,  (1  •num_src) ) ' ) )  *ccnj  (DS) ))]'). \  (N*N*anes  (M,  1) ) )  ; 

%  NORMALIZE  RESPONSE,  CONVERT  TO  DB 

br_plot  ( : ,  k)  =  10 .  *logl0  (br_plot  ( : ,  k) .  /max.  (br_plot  ( : ,  k) ) )  ; 

%  ADD  FACTOR  TO  EACH  RESPONSE  FOR  DISPLAY  PURPOSES 

br_plot(:,  k)  =  br_plot(:,  k)  +  (k*20)  ; 
and; 

%  PLOT  MUSIC  BEARING  RESPONSES 
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plot  (bearing_angles ,  br_plot ,  '  cl  -  ’ )  ; 


%  Label  x  axis 
%  label  y  axis 


xlabel ( ' cos (bearing) ' ) ; 
ylabel ( ' Iteration  ; 

%  ******************************************************** 
%  FUNCTIONS  -  MUST  BE  UTILIZED  AS  SEPARATE  M-FILES 
%  ******************************************************** 

function  y  =  A1_F (a,  N)  ; 

%  COMPUTES  THE  MATRIX  Al 

y  =  toeplitz(  [l;a;zeros(  (N-length(a)-l) ,  1)],  ... 

[1  zeros  (1,  N-l)  ] ) ; 

return; 

%  ******************************************************** 

function  y  =  A3_F(a,  N) ; 

%  COMPUTES  THE  MATRIX  A3 

y  =  toeplitz  ( [zeros  (N-length(a) ,  1) ;  flipud(a)  ] ,  zeros  (1,  N) ) ; 
return; 

%  ******************************************************** 

function  y  =  ARMMa,  b,  theta,  sp)  ; 

%  COMPUTES  SPATIAL  DENSITY  OF  AN  ARMA(P,  Q)  PROCESS  ACROSS 
%  BEARINGS  SPECIFIED  BY  THETA,  AND  AT  SENSOR  SPACING  SP 

p  =  length  (a) ;  %  AR  order 

q  =  length(b) ;  %  MA  order 

M  =  length  (theta) ;  %  Manber  of  bearing  bins 

z  =  exp(-sqrt  (-1)  *2*pi*sp*cos(theta) )  ; 

A  =  a*  [  [ones  (p,  l)*z]  .A[  [  [0:1;  (p-1)  ] .  ’  ]*ones(l,  M)  ]  ] ;  %A(z) 

B  =  b* [  [ones (q,  1)  *z]  .A[  [  [0:1:  (q-1)  ] . '  ]*ones(l,  M)  ]  ] ;  %  B(z) 
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y  =  (B.*conj  (B) ) .  /  (A.  *conj  (A) )  ; 


%  Spatial  density 


return; 

%  *********************************************** 

function  y  =  B_F(vamoiz/  a,  N)  ; 

%  COMPUTES  COVARIANCE  MATRIX  BASED  ON  THE  GOHBERG 
%  FORMULATION 

A1  =  Al_F(a,  N)  ; 

A3  =  A3_F(a,  N) ; 

%  True  noise  covariance  matrix 

y  =  inv(  (1/vamoiz)  *  (Al*Al '  -  A3* A3 ' ) )  ; 

return; 

%  ******************************************************** 

Function  y  =  B_FI(vamoiz,  a,  N) ; 

%  COMPUTES  INVERSE  COVARIANCE  MATRIX  BASED  ON  THE 
%  GOHBERG  FORMULATION 


%  Oczqpute  matrix  A1 
%  Ccmpute  matrix  A3 


Al  =  Al_F(a,  N) ;  %  Caqpute  A1  matrix 

A3  =  A3_F  (a,  N) ;  %  Ocnpute  A3  matrix 

%  Inverse  of  true  noise  covariance  matrix 

y  =  (Al*Al '  -  A3* A3 ' )  /vamoiz; 

return; 

%  ******************************************************** 

function  y  =  Z_FUNC(i,  N)  ; 

%  COMPUTES  THE  Z(I)  MATRIX 

y  =  zeros  (N)  ; 

jjnatrix  =  (1:N) . ’♦anesd/N) ; 
this  =  find(  (j_matrix-j_matrix. ' )  ==  i); 
y(this)  =  onesd,  length  (this)  ); 
return; 
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